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Abstract. Theoretical understanding of strongly correlated systems in one spatial dimension (ID) 
has been greatly advanced by the density-matrix renormalization group (DMRG) algorithm, which is a 
variational approach using a class of entanglement-restricted states called Matrix Product States (MPSs). 
However, DRMG suffers from inherent accuracy restrictions when multiple states are involved due to 
multi-state targeting and also the approximate representation of the Hamiltonian and other operators. 
By formulating the variational approach of DMRG explicitly for MPSs one can avoid errors inherent 
in the multi-state targeting approach. Furthermore, by using the Matrix Product Operator (MPO) 
formalism, one can exactly represent the Hamiltonian and other operators relevant for the calculation. 
The MPO approach allows ID Hamiltonians to be templated using a small set of finite state automaton 
rules without reference to the particular microscopic degrees of freedom. We present two algorithms 
which take advantage of these properties: eMPS to find excited states of ID Hamiltonians and tMPS for 
the time evolution of a generic time-dependent ID Hamiltonian. We properly account for time-ordering 
of the propagator such that the error does not depend on the rate of change of the Hamiltonian. Our 
algorithms use only the MPO form of the Hamiltonian, and so are applicable to microscopic degrees of 
freedom of any variety, and do not require Hamiltonian-specialized implementation. We benchmark our 
algorithms with a case study of the Ising model, where the critical point is located using entanglement 
measures. We then study the dynamics of this model under a time-dependent quench of the transverse 
field through the critical point. Finally, we present studies of a dipolar, or long-range Ising model, again 
using entanglement measures to find the critical point and study the dynamics of a time-dependent 
quench through the critical point. 



1. Introduction 

The great success of experimental ultracold atomic physics has made the study of strongly correlated 
one-dimensional (ID) quantum systems a major avenue of current physics research. Examples of 
^ ■ novel condensed matter physics realized with ID atomic systems include the role of integrability in 
thermalization pQ and static [2] and dynamic j3] quantum simulators of Hubbard models. Furthermore, 
as ultracold molecules approach quantum degeneracy [H El E], lattice models with complex internal 
degrees of freedom and long-range interactions become relevant El E]- ^ s more an d more complex 
models become amenable to study, the need for numerical methods which can adapt to different degrees 
of freedom, different Hamiltonians, and different dynamical processes thus becomes essential. 

In addition to practical interest in understanding and benchmarking atomic and molecular quantum 
simulators, the ability to simulate the dynamics of ID systems also provides insight into fundamental 
questions such as the universality of dynamics approach quantum critical points and the effects of 
integrability on the thermalization process [10] . The natural setting for studying dynamics near critical 
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points is a quantum quench where one of the parameters of the Hamiltonian is driven through a quantum 
critical point following a time dependent protocol, for example 

g(t) =9o+ v{t ~ to)r o(t-t ) , (i) 

with 6 {t) the step function. Such quenches pose a special difficulty for numerical studies as by definition 
they involve evolution with a time-dependent Hamiltonian which does not commute with itself at 
different times. The propagator is then generally a time-ordered exponential whose precise form may 
be difficult to ascertain. Standard methods such as the Suzuki- Trotter expansion which ignore the time 
dependence of the Hamiltonian require [IT] that the infinitesimal time step used be much less than 
the fluctuation time-scale of H (t) to be valid, 5t <C \dH / dt\~ . This can cause simulations with rapid 
quench rates to become numerically very costly, and invalidates the approach altogether for non-analytic 
time dependence. 

Currently the only unbiased method available for the dynamics of quantum systems is exact 
diagonalization (ED). By unbiased, we refer to the fact that the other methods available for dynamics 
are generally variational, and so have a bias towards a particular ansatz. ED is limited in an essential 
way by the exponential growth of the size of the Hilbert space with the physical size of the system. 
The current state of the art is ~ 40 spins for spin- 1/2 models and 20 sites for a fermionic Hubbard 
model at half filling. These sizes are often too small for accurate finite-size scaling. An extremely 
powerful method for the low-energy properties of ID systems is White's Density Matrix Renormalization 
Group (DMRG) algorithm, which uses a variational ansatz based on a class of states known as Matrix 
Product States (MPSs). MPSs will be reviewed in Sec. 12.11 DMRG uses an iterative procedure to 
develop a set of reduced bases that the full many-body problem is projected into, and then variationally 
minimizes an energy functional in this reduced space, enlarging it if necessary. DMRG uses an implicit 
MPS representation, which is to say that the state is not stored explicitly. This also means that the 
Hamiltonian and other operators in the calculation are stored in an approximate way, as they are 
represented within the reduced basis describing the variational state. This does not cause problems in 
practice when a single state is sought using the DMRG process. In fact, one can show that the algorithm 
to variationally find the ground state is identical when phrased in the implicit formulation of DMRG 
and when using an explicit MPS representation for the variational state, other than the representation 
of the Hamiltonian [12]. However, because of the exact representation of operators independent of the 
state, MPSs can put rigorous bounds on distances from exact quantum states such as eigenstates by 

considering quantities like the variance — E^j \ip) with H the Hamiltonian operator and E the 

energy expectation of the MPS In contrast, DMRG can only return the distance of the variational 
state from the approximate operator H in the given variational basis, and is unable to determine how 
well the given variational basis represents the true operator. A particularly clear indication of the 
failures this can cause is given in Ref. (13] where time evolution of a particular initial state in DMRG 
fails because the Hamiltonian has no nonzero matrix elements in the initial DMRG basis. 

The situation becomes much different when multiple states are sought using the DMRG procedure. 
In this case the reduced density matrix used to determine the optimal reduced bases for the algorithm is 
a convex sum of the reduced density matrices for the desired states. This is called multi-state targeting. 
In contrast, an explicit MPS representation stores each of the desired states separately as an MPS. In 
multi-state targeting, none of the states can be represented with the same accuracy available if DMRG 
targeted that state alone. The MPS representation also deals automatically with the fact that each 
state has its own optimal bases for representation, whereas in DMRG these bases are all tied together 
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by the multi-state targeting. In this work we present two algorithms which take advantage of MPSs' 
ability to deal with multiple states, eMPS to find excited states of ID Hamiltonians and tMPS to 
simulate the dynamics of a generic time-dependent Hamiltonian. In the first algorithm a projector 
orthogonal to a set of lower-lying eigenstates is constructed from their MPS representations and used to 
orthogonalize a variational state against this set. In the second algorithm Krylov vectors in a Lanczos 
approximation to the matrix exponential are stored separately as MPSs and combined in an optimal 
way only at the end of the calculation. While Krylov-based MPS approaches have been used [331 166] 
to study time-dependent processes, the errors in these approaches were set by time derivative of the 
Hamiltonian. In contrast, by taking explicit account of the time ordering of the propagator, the errors 
in our approach are set only by commutators of the Hamiltonian at different times, and hence allow for 
larger time steps. Because of the explicit MPS representation, we are able to put bounds on the errors 
in each step of the calculations. 

Finally, MPSs have a natural operator-valued extension known as Matrix Product Operators 
(MPOs) which allow for the exact representation of all operators used in the calculation. We present 
a general framework for constructing MPOs from a set of rules which is independent of the nature 
of the microscopic degrees of freedom. This allows for the templating of ID Hamiltonians for general 
purpose software. In addition, the ability to perform arithmetic operations on MPOs exactly enables 
us to perform time-evolution using our tMPS algorithm with knowledge only of the MPO form of 
the Hamiltonian and the time-dependent functional form of its parameters. To emphasize the general 
nature of our algorithms, we include a generic simulation protocol for the out-of-equilibrium dynamics 
of strongly correlated ID systems using the algorithms presented in this paper. 

The remainder of this paper is organized as follows. In Sec. [2] we review the theory of MPSs, 
MPOs, and their canonical forms. In addition to providing a canonical form for operators within the 
matrix product formalism, we define finite state automaton rules for MPOs and demonstrate how ID 
Hamiltonians can be constructed from a small set of such rules. In Sec. [3] we review the algorithm 
for finding ground states using MPSs as variational ansiitze, and in Sec. H] we present the eMPS 
algorithm which extends the ground state search to general excited states. Sec. [5] discusses how to 
extract observable quantities from MPSs. In Sec. |6]we discuss methods for time evolution with MPS. In 
particular, we provide the tMPS algorithm to time evolve an MPS using only the MPO representation 
of a Hamiltonian and the functional form of its time- dependent parameters. We contrast our approach 
with other Krylov subspace approaches and identify the possible sources of error. In Sec. [7] we present 
two case studies. The first is of the Ising model in a transverse field, where we study both the statics and 
the dynamics of a linear quench of the transverse field through the quantum critical point. The second 
is of a dipolar, or long-range Ising model in a transverse field, where we also determine the critical 
point from the statics and study a linear quench of the transverse field. Finally, in Sec. El we conclude. 
Details concerning numerically exact solutions for the Ising model which are used to benchmark our 



algorithms are given in |Appendix A 



2. Brief review of Matrix Product Formalism 



2.1. Matrix Product States 

The Hilbert space of a quantum mechanical many-body system is exponentially large in the physical 
size of the system, for example the number of unit cells in a lattice or the number of particles. Stated 
another way, a state picked at random from the Hilbert space of a quantum many-body system will have 
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entanglement (as quantified by the Schmidt measure [H]) which grows exponentially with the system 
size. In contrast to this random state, it has been shown that the class of states which are physically 
relevant in the sense that they can be prepared from some reference state by generic time evolution in 
polynomial time [IT] or are useful for quantum computation [T5] is much smaller than the full Hilbert 
space. In ID, the physically relevant class of states appears to be those which have entanglement which 
is either constant or polynomially growing as a function of system size. A convenient representation of 
states with entanglement restricted in this manner is known as matrix product states (MPSs) [T6l IT71 IT2] . 

We consider our physical system to comprised of a ID lattice of sites, where each site i is a d- 
dimensional Hilbert space Hj spanned by the vectors We will refer to d as the local dimension, 

and take all sites to be identical for simplicity. We define an MPS on a lattice with L sites as 

d 

|^mps)= Tr{A^...A^)\h...i L ), (2) 

h,...i L =l 

where the object A™ lk is a Xk x Xfc+i matrix (with xi — Xl+i) and Tr denotes the matrix trace. We will 
refer to the maximum linear dimension of any of the matrices A^ lk , max^Xk, as the bond dimension 
of the MPS, and denote this quantity by x- X niay be used as a measure of the entanglement of the 
state [12]. In this work we will focus on systems with open boundary conditions (OBC). MPS algorithms 
can also be devised for systems with periodic boundary conditions, as discussed in Refs. [HI [191 13 I2H] , 
but these algorithms have worse scaling and are generally less numerically stable than their OBC 
counterparts. For OBC, Xi — Xl+i = 1, and arguments using the Schmidt decomposition demonstrate 
that Xk < mm(d k -\d L - k ) [2TJ. 

MPSs have been used for many years to represent exact ground states of parent Hamiltonians [22] 
which are formed from projectors onto local high-symmetry subspaces [231 121] • However, it was not 
until the great success of the density matrix renormalization group algorithm (DMRG) pioneered by 
White [22] that MPSs became valuable as variational ansdtze in their own right [18]. Why are MPSs 
useful as variational ansatzel It has been shown [26] that the ground states of gapped ID systems have 
bipartite entanglement which does not depend on the system size. Such states can be represented exactly 
as MPSs with a fixed bond dimension [27]. This is an example of an area law [21]; the entanglement 
between two disjoint subsystems depends only on the boundary of the two regions and not on their 
volume. For systems near a quantum critical point which is described by a conformal field theory 
(CFT), this area law is subject to weak logarithmic violations such that the entropy of entanglement 

between two subsystems of size L is given by the Calabrese-Cardy formula [211 130] 

c 

S L ~ a + - logL, (3) 
o 

where a is a constant and c is the central charge of the underlying CFT. Here ~ denotes scaling 
equivalence in the bulk of an infinite system. In finite systems there are often oscillating boundary 
and finite size contributions [311 E2]- Hence, the bond dimension of an MPS describing a conformally 
invariant critical system is given as xl ~ exp S l ~ e a L c / 6 , which grows only polynomially in the system 
size. Typical values of c range from 1/2 for the Ising model [33] and 1 for the Bose-Hubbard model [M] 
to 2 for more exotic phases like the gapless Mott insulator of the JK model [35]. We note that, strictly 
speaking, finding an MPS which approximates the ground state of an arbitrary ID Hamiltonian to 
an accuracy which is an inverse polynomial in the system size is still NP-complete [36], but practical 
experience demonstrates that this method is extremely useful and robust for physical systems of interest. 

We adopt the following conventions for the representation of tensors: we use roman indices for 
indices which correspond to physical states and greek indices for indices which are summed over in the 
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(a) (6) (c) 



Figure 1. Diagrammatic notation for tensor networks, (a) A rank-three tensor is represented as a 
point with three lines extending from it. (b) Contraction of two rank-three tensors is accomplished by 
connecting the contracted index, and produces a tensor of rank four, (c) An MPS on 5 sites with open 
boundary conditions 

is represented as a contraction over rank three tensors with two rank-two boundary tensors. 



matrix-product ansatz. Explicitly writing out Eq. ([2]) with these indices, we have 

I^mps)= £ £ A^...A^_\ ai \i x ...i L ). (4) 

oi-ai-i u,...il=1 

A superscript index in square brackets [] denotes the lattice site that the physical indices of the tensor 
describe. A superscript index in curly braces {} denotes association with a particular many-body state. 
For example, the MPS tensors at site j of the MPSs \4>k), k = l,...,n would be denoted A^ k ^ J . 
Finally, indices which appear together in parentheses, e.g. (a/3), represent a composite index which 
runs over the Cartesian product of the indices in the parentheses. As an example, if a = 1, . . . , Xa and 
j3 = 1, . . . xp> ( a fl) = (a — l)xp + P- To lighten the notation, we will leave off indices when they are 
unnecessary. 

A particularly useful means to visualize MPSs and manipulations with them is provided by tensor 
network diagrams like those shown in Fig. [TJ [37]. Here a rank-/c tensor is represented by a point with k 
lines extending from it. Each line represents one of the indices of the tensor. Whenever a line connects 
two points, that index is summed over, and disconnected lines represent free indices. Hence, an MPS 
can be represented as a chain of rank-3 tensors as in Fig. QTc). Note that the first and last MPS tensors 
are rank two because we have assumed OBC and so xi — Xl+i = 1- 

We note that the MPS definition Eq. (J2J) does not uniquely specify the tensors A. That is, we can 
insert an invertible matrix X and its inverse A -1 between any two MPS matrices without altering the 
physical content of the state: 

d 

|^mps)= Yl Tr(AWi...AW*)\i 1 ...i L ), (5) 

h,...i L =l 

J[\j]ij = aWjx , A^'+ife+i = x~ 1 A^ j+1 ^ i+1 , (6) 

d 

I^mps) = Tr { All]h ■ ■ ■ A^A^ +1 ^ . . . A^) i L ) = I^mps) • (7) 

ii,...i L =l 

This is referred to as gauge freedom in the literature [21]. For OBC, we can specify the state uniqueljlj] 
by choosing a site k, which we call the orthogonality center of the MPS, and requiring that all sites % 



% The state is unique up to possible degeneracies in the Schmidt decomposition. 



OOE dynamics with MPS 



6 



(«) (b) ( C ) 

EC=C =m 

Figure 2. Gauge conditions for MPSs. (a) Left gauge condition Eq ©. (b) Right gauge condition 
Eq ©. (c) Using (a) and (b), the norm squared of an MPS reduces to the trace of its orthogonality 
center squared, here chosen to be the third site. 



to the left and right of k satisfy the left 

^2A i] A l = I (8) 

and right 

A'A i] = I (9) 

i 

gauge conditions, respectively. In these expressions, / is the appropriately dimensioned identity matrix. 
These conditions are shown in graphical notation in Fig. [21(a) and (b), respectivelyl§] Graphically it is 
clear that the norm squared of the state is 

(^) = ^Tr (A^A^ k A , (10) 

i 

as shown in Fig. [2(c), and so this site carries all information about the norm of the state. This particular 
canonical form for an MPS is called mixed canonical form [12]. The mixed canonical form is crucial for 
improving the speed and numerical stability of variational algorithms with MPSs. 
We can impose the left gauge conditions via the following recursion: 

AiaQfi =48, (11) 

^ f/ (ai ) 7 S 7 V r 7 t /3 = A^p , (12) 

7 

A [ S =U {mh , (13) 

=E s ^VC^ (14) 

p 

where Eq. ffl2|) represents the singular value decomposition (SVD) of A with E a diagonal matrix of 
singular values and U and V unitary matrices!]]] Because U returned from the SVD is unitary, Eq. (JSJ) 
is satisfied by construction. Similarly, the recursion for the right gauge conditions is 

Mm = a { % , (15) 

u <*i^i v kp) = Mw ' ( 16 ) 

7 

a 

§ Here we also establish the graphical convention that downwards pointing lines correspond to Hermitian conjugates of 
tensors. 

|| Note that any matrix decomposition of A which returns a unitary matrix as part of the decomposition will suffice in 
place of the SVD. In particular, the QR decomposition [35] is particularly efficient when the rank of A is not required. 
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(a) 



Figure 3. MPO in diagrammatic notation, (a) A matrix product operator consists of a contraction of 
rank-four tensors. We adopt the graphical convention that the line below the horizontal corresponds to 
%' and the line above the horizontal to i for an operator Ou> which takes the local state from \i') to 
(b) The product of an MPO and an MPS produces another MPS whose bond dimension is the product 
of the bond dimensions of the original MPS and the MPO. 



To put a general state into mixed canonical form with orthogonality center k we begin at site 1 and 
iterate Eqs. (|Tl| -( |T4|) until we reach site k, then start at site L and iterate Eqs. ( TT5|) -( TT8|) until we again 
reach site k. 

Finally, we note that the set of all MPSs with a fixed bond dimension \ is not a vector space, as 
the sum of two MPSs with bond dimensions xa and xb has a bond dimension x which is bounded by 
the sum of the two bond dimensions x < Xa + Xb- This can be seen from considering the sum of the 
two states |0 . . . 0) and |1 . . . 1), with MPS representations 

0...0) Y ':•••'/■>• ^ =(£,„), (19) 

h...i L 

|l...l) = Y B*...B i *>\i 1 ...i L ), B i = (5 i: i) . (20) 

The matrices A 1 and B l have a bond dimension of 1, as these are product states. Their sum is 

|0...0) + |l...l) = Y C [1]il C l2]i2 . . .C [L - 1]iL ~ 1 C [Lfe |z 1 . . .i L ) , (21) 

h—iL 

CW=(5 i0 S a ); CW=( S £ ^J,2<j<L-l; C™ = ^ ); (22) 
which has a bond dimension of 2. 



2.2. Matrix Product Operators 

The natural operator generalization of MPSs is a Matrix Product Operator (MPO), defined as 

d d 

Ompo= Y E (23) 

il,...i£=l i' ly ..i' L =l 

Here W^ u ' is a Xk x Xk+i dimensional matrix, and we will again refer to the maximum value of x° as 
the bond dimension of the operator. Note that this bond dimension x° need not be the same as the 
bond dimension x appearing in the MPS representation, Eq. (J2J). That an MPO takes MPSs to MPSs 
can be seen clearly from the graphical representation of Fig. [3j We also see from this representation 
that the bond dimension of the MPS representing the product of an MPO and an MPS is the product 



OOE dynamics with MPS 



8 



of the bond dimensions of the MPO and the MPS. Because the states \i\ . . . il) are tensor products, we 
can also use the notation 

d d 

MPO = Tr . . . W [L]lLl '^ , (24) 

h,...i L =l i' v ...i' L =l 

where each one of the objects W [i M = W^'^ij)^. That is to say, we can consider the matrices 
which appear in the matrix-product ansatz of an MPO to be operator-valued. 

It is tempting to look for canonical forms for MPOs just as we did for MPSs, but the relevant 
norm for MPOs is the Frobenius norm (Oi,02) = r Ti{0\C>2) which scales exponentially in the local 
dimension with the number of lattice sites. Thus, for a typical many-body system with an exponentially 
large Hilbert space, the elements of the orthogonality center can differ in magnitude greatly, causing a 
catastrophic loss of precision during orthogonalization. However, most physically relevant MPOs such 
as one- dimensional Hamiltonians can be written down exactly in terms of an MPO canonical form which 
is analogous to an LU decomposition [39]. To discuss this canonical form, it is useful to recast an MPO 
as a finite state automaton (FSA) [4"0| |4T]. 

To recast an MPO as a FSA, we first enumerate all of the physical operators we use to define our 
local Hilbert space, O = {O a }, where a labels distinct operators. We will call this set our operator 
alphabet (OA). As examples, the OA for the Ising model would be j J, a x , <j 2 j with J the identity 
operator and a v the Pauli operator along the u th Cartesian direction. The OA for the Bose-Hubbard 
model (121 H3] would be j/, & , b, nj, where b is a bosonic destruction operator and n = &b the number 
operator. The particular matrix representation of the OA fixes the local basis states {\i)}. 

Using the OA, we now introduce a set of FSA rules {7Z P (^{0 Pl , . . . , Pn }, {h p } , w p ^j } which generate 

strings of the n operators {0 Pl , . . . , Pn } G O weighted by scalar w p and variables {hp}. We will call 
the variables {h p } Hamiltonian parameters and the scalars w p weights. Each rule consists of three 
operator- valued matrices in the case where the {h p } do not depend on position. The first matrix is the 
rightmost matrix in the MPO representation, W^ L \ and represents the initial configuration. Next, we 
have the MPO matrix of the bulk of the chain W^, 2 < j < L — 1, which takes an input vector of 
operators on k sites and produces a vector of operators on k + 1 sites according to some deterministic 
pattern. Finally, we have the leftmost matrix in the MPO representation, W^, which extends the 
operators according to the patten of and returns a 1 x 1 operator-valued matrix. The trace of 
this matrix as in Eq ( 123]) is the desired Hamiltonian term. The generalization to position-dependent 
Hamiltonian parameters requires L — 2 matrices in place of the bulk matrix W^, 2 < j < L — 1, but 
the only modification is that h p becomes h p (j). 

As a concrete example, consider the Ising model 

H=-J^2af^-h^2af. (25) 

The Hamiltonian consists of two rules. The first is a site rule TZ S i te (a x , h,—l) which generates the string 
—hY2i&f. The three matrices which provide this rule are 

nfi = ( / ) , W JL— = ( _i, ° ) , nfil = ( _l, ) • («) 



OOE dynamics with MPS 
As can be verified, 

k 

n wsxs 



j=L-l 



I. ..I 



(27) 



and so this rule produces the desired operator. Similarly, there is a bond rule TZbondiW 2 , <? z }, J, — 1) 
given by 



bond 



( i o o \ 

<T 2 

-Jb z I 



bond 



/ J \ 

V / 



(2« 



which produces — JJ2(ij)^i^j^ with denoting a sum over nearest neighbors i and j. The full 

Hamiltonian is given by the direct sum of the matrices. Collecting rows of the direct sum which are 
exactly parallel, we have the MPO representation of the full operator 

/ / \ 



'^'ising 



-ha x -Ja z I 



[2<j<L-l] 
Ising 



V 



a z 
-ha x -Ja z I 



'^'ising 
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( 1 \ 

a z 

l — ha x j 



(29) 



This construction can be readily extended to general sums of nearest-neighbor and on-site interactions 
with Hamiltonian 

n b n s 

h = - E J « E d ^ d f a - E h p E 6 ? ( 30 ) 

a=l 0=1 i 

by summing up the individual rules: 



yy[2<j<L-i] 



J 

6 62 -i> o 



o 
o 



o\ 





-Jn b O b ^ I 



w [L] 



J 



I 

d b2 - 



\ 



(31) 



.(32) 



The bond dimension of the MPO representation of this Hamiltonian is 2 + ub- MPO representations 
are not restricted to nearest-neighbor terms. Exponentially decaying terms of the formal 



H = -jJ2^ XU ~ l ' 1] OlO] 

i<j 

can also be accommodated with the rule 7^. e xp({0 1 , O 2 }, {J, A}, — 1) 



(33) 



exp 







/ 


i 








\, 




/ 




\ 


o -jo 1 i ) 


W [2<j<L-l] _ 
' ' K exp 




6 2 









w lL] = 

' ' K exp 




6 2 






V 





-JO 1 


J 


J 


V 





/ 



We can also produce a general monotonically decaying term 
H = jJ2fti 



- i) 0\0] 



(34) 



(35) 



% Note in this form that the nearest neighbor coupling is J and only longer ranged couplings contain A. 
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by approximating the term as a sum of exponentials to a desired distance r cuto g. That is, we minimize 
the functional 



'exp ' cutoff' 



F(a,b) = ^ Yl \f( x )~ a i b i 



-H 2 



(36) 



with respect to the exponential weights a« and decay parameters 6j, where a and b are the elements 
Oj and bi, respectively, arranged as vectors and n exp is a convergence parameter controlling the number 
of exponentials used in the expansion [HJ H3 HSj. Note that the fit is only guaranteed to be 
accurate to r cuto a while the term has infinite range. This does not cause difficulties in practice for 
decaying functions, as the resulting fit is also guaranteed to be monotonically decaying. The decaying 
function rule is then just a sum of these n exp exponential rules, TZdi ({O 1 , O 2 }, {J, f(x),n exp , r cuto fj}, 1) = 
Er=T ^({O 1 , O 2 }, {Ja h - log&i}, 1). In contrast to the other rules presented above this rule is 
not exact, but the number of exponentials can be increased to any desired accuracjQ-a few tens of 
exponentials suffice for a 1/r 3 interaction on a lattice of 50 sites [35] . While the bond dimension 
increases linearly with the number of exponentials, the number of nonzero terms in the MPO also 
grows only linearly with the number of exponentials, and so operations with the MPO scale well as n cxp 
increases. However, the eigenstate of a Hamiltonian with larger n exp may be more highly entangled 
than with smaller n exp , leading to an increase in x and longer runtimes. That is to say, the relationship 
between the bond dimension of an MPO and the bond dimension of an MPS representing an eigenstate 
of the MPO is difficult to predict. 

While the pure functional interaction Eq. ( |35l) is appealing from a theoretical point of view, in 
practical applications infinite range interactions do not represent a consistent level of approximation. 
That is to say, at some distance interactions are screened, where the screening length is set by energetics 
or lifetime constraints in the case of ultracold atoms or molecules [U]. Hence, we provide the finite- 
ranged rule T^FimteFunctiont^O 1 , O 2 } , {h, f, r cutoff }, w) which generates the Hamiltonian 



j+^cutoff 



j=i+l 



and is given by 

= ( o o ... o ... o whf(i)6 l i 
( i o ... o 



(37) 



(31 



yy[2<j<L-l] 



0> 



o 2 





/Qcutoff) 
/Ocutoff-1) ' 



6 2 





o 
o 

m 
f(k-i) 















o\ 






w [L] 



( i \ 

6 2 
6 2 

6 2 



6 2 
o 



• (39) 



V / 



... ... jffil 

y o o ... o ... o whf(i)6 l i J 

For a vector / which contains zero elements a small modification must be made, but the bond dimension 
remains the same. 

A small set of rules like those presented above allow us to template operators for MPS simulations. 
That is, given a method to produce a Hamiltonian from a set of basic rules, we can generate a wide 

+ By accuracy we mean that the functional Eq. (|36[) is smaller than a given tolerance. 
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variety of different Hamiltonians by specifying different OAs, weights, and Hamiltonian parameters 
as inputs to a program. The resulting MPOs are lower triangular and usually very sparse, and so 
sparse matrix structures can be used for efficiency. Beyond the flexibility that MPOs provide to MPS 
algorithms, they also allow for arithmetic operations such as addition and multiplication to be performed 
exactly, albeit at the expense of a growing bond dimension. The sum of two operators expressed as 
MPOs has MPO matrices which are the direct sums of the constituent matrices, and the product of 
two MPOs has MPO matrices which are the direct product of the constituent MPO matrices. Both of 
these operations preserve the sparse lower triangular structure of the MPOs. 



3. Variational Ground State Search 



We now turn to using MPSs as variational ansatze for the eigenstates of a Hamiltonian expressed as an 
MPO. The ground state is found by minimizing the functional 

£U)} = (^\H\^)-E(^) (40) 

with respect to the parameters of where E is a Lagrange multiplier enforcing normalization. The 
general minimization of this functional is an NP-hard problem, so we instead adopt a local search 
heuristic that has proven to work well in practice. Let us consider fixing all parameters in the MPS 
except for a contiguous block of s MPS tensors A® . . . A^ +s ^. The single-site (s = 1) and two-site 
(s = 2) algorithms are the most commonly used variants. We then find the extremum as 

which corresponds to the diagrammatic equation shown in Fig. HI Here, the partial derivative with 
respect to a tensor is defined to be a tensor whose elements are the partial derivatives with respect 
to the elements of the tensor. If we assume that the block of tensors to be optimized contains the 
orthogonality center then the rightmost diagram reduces to the block of tensors being optimized and 
the leftmost diagram is the action of the effective Hamiltonian on this same block. Thus, minimization 
consists of finding the eigenvector corresponding to the smallest eigenvalue of the effective Hamiltonian 
eigenvalue problem 

H^A® . . . A j+S ~^ = EA® . . . A^ +s ~^ , (42) 
where the effective Hamiltonian is 

* Wi-h+.-ii'j-tj+s-l _ ST V 7-0 W Wii * H7 b'+s-l]ij+»-iij+,-in[,-+ S ] (A o\ 

kk' 71... 7s 

L and R are the partial overlaps of the Hamiltonian MPO with the state as in Fig. H and is the 
MPO tensor at site j of the Hamiltonian. represents the Hamiltonian for the variational degrees 
of freedom in the block to be optimized with the rest of the state held fixed. This justifies our use of E 
as the eigenvalue, as E obtained from the solution of this equation is the current best estimate for the 
energy. We can view Eq. ( 1421) as a linear eigenvalue problem by combining indices as 



Yl ^cff {ai J ...i J+a _ 1 p)( a >i' j ...i' j+B _ 1 P>) ■ ■ ■ Ab+S ^ (a'i' r ..i' j+s _^') 

The linear dimension of this matrix representation of H^q is Xjd s Xj+s-ii an d so a solution of this problem 
with dense methods would require 0(x^d 3s x^+ s -i) basic operations, leading to a very slow algorithm 
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Figure 4. Variational ground state search in diagrammatic notation. Here we display the single site 
(s = 1) effective Hamiltonian eigenvalue problem for simplicity. The contractions between the block of 
tensors to be optimized, the MPO, and the rest of the diagram have been omitted to accentuate the 
structure of the effective Hamiltonian. The rightmost equality follows from assuming that the variational 
site is the orthogonality center. 



of order 0(x 6 )- In contrast, by taking advantage of the separable form of the effective Hamiltonian 
Eq. (T4"3"j) multiplication of our block of tensors by the effective Hamiltonian can be done in 0(x 3 ) 
time [12] . Thus, sparse eigensolvers such as the Lanczos [38] or Davidson [38] algorithms, which require 
only matrix-vector multiplies, should be employed to solve this eigenvalue problem. 

The general algorithm for ground state search is thus as follows. We begin with an initial state with 
orthogonality center at site k. We choose a block of tensors containing k and optimize them by solving 
the effective Hamiltonian eigenvalue problem. We then shift the orthogonality center and the block of 
tensors one site to the right and again optimize. We continue shifting to the right until we reach the 
right boundary. We then reverse direction, shifting the orthogonality center and the block of tensors to 
be optimized to the left and solving the effective Hamiltonian eigenvalue problem until we reach the left 
boundary, at which point we reverse direction again. A single iteration of this optimization cycle which 
affects each tensor twice is called a sweep, and sweeping is continued until convergence. In addition, 
using the MPO form of the Hamiltonian, it is possible to develop a caching algorithm for the overlaps L 
and R such that the solution of this problem requires O (L) scaling in the number of lattice sites |41jl*l 
To do so, we begin the iteration with a guess for the ground state assumed to have orthogonality 
center k. We then use the left recursion 



■ [m] 



E 



E 



k'k 



aW 



to generate the L overlaps up to k and the right recursion 

/ j kk: 



R 



k/3'P 



2^ 



V R [e+1] A 



(45) 
(46) 

(47) 
(48) 



to generate the R overlaps down to k + s. Here the square braces indicate the order in which the 
contraction should be performed to achieve ideal scaling. Once the eigenvalue problem has been solved 
and the orthogonality center shifted, we use the recurrence Eq. PSI) to update the overlaps when we 



* This scaling does not account for possible L dependence of the bond dimension \ such as exists for conformal critical 
points. 
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are sweeping to the right and the recurrence Eq. (148]) to update the overlaps when we are sweeping to 
the left. 

Convergence is achieved when the variance 

A = (^-E 2 )), (49) 

with E the energy eigenvalue, drops below a user-specified tolerance e. Given the MPO form of the 
Hamiltonian, the variance operator A = H 2 — E 2 can be constructed by constructing an MPO whose 
matrices consist of the direct product of the corresponding matrices from H, = ® 
and then subtracting —IE 2 /L from the lower leftmost element of each where / is the identity 

operator. This representation is exact, in contrast to DMRG-based approaches where the basis of the 
Hamiltonian is tied together with the basis of the state. The variance is a much better measure of 
convergence of the state than the so-called discarded weight which is used to measure convergence of 
the two-site DMRG algorithm. This is because it is a property of the actual MPS state and not of the 
eigenvalue. As a note of caution, the variance only guarantees that the state found is an eigenstate 
to the given tolerance, it does not specify that it is the ground state. This has not proven to cause 
problems in practice for non-disordered systems. 

In summary, the complete algorithm for variational ground state search is: 

(i) Input: Input a Hamiltonian H in MPO form, an initial guess for the ground state in MPS 
form, and a tolerance e for the variance. 

(ii) Initialization: Construct the LR overlaps using the recursions Eq. ( H6|) and ( 148]) . 

(iii) Sweeping : 

(a) Solve the effective Hamiltonian eigenvalue problem Eq. ( |42|) and replace the variational block 
of tensors with the eigenvector corresponding to the lowest eigenvalue. 

(b) Shift the orthogonality center to the right and update the LR overlaps using the recursion 
Eq. flUD. 

(c) Continue to iterate a and b {right sweeping) until the right boundary is reached. 

(d) Solve the effective Hamiltonian eigenvalue problem Eq. ( 142|) and replace the variational block 
of tensors with the eigenvector corresponding to the lowest eigenvalue. 

(e) Shift the orthogonality center to the left and update the LR overlaps using the recursion 
Eq. (08]). 

(f) Continue to iterate d and e {left sweeping) until the left boundary is reached. 

(g) Iterate a and b until k is reached. 

(iv) Check convergence: Using the most recent estimate of the energy eigenvalue E from the last 
effective Hamiltonian solution, construct the variance operator A and find the variance. If A < e, 
exit, otherwise return to a. 

4. Variational Excited State Search 

We now turn to finding excited states. We find the n th excited state by minimizing the functional 

n-l 

£ [|V)] = (V|%> - £(V#> - J>*W*> (50) 

k=0 

where {|0fe}} are the n — 1 lower-lying eigenstates of H and the {/!&} are Lagrange multipliers enforcing 
the orthogonality constraints (ip\<ftk) = 0. Again fixing a block of s tensors, the minimization of this 
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£[j]{k} ftti+sHk} 



Figure 5. Linear forms enforcing orthogonality, (a) The overlap (tp\<fik) m diagrammatic notation. 
The thin red lines correspond to (ip\ and the thick black lines to \<f>k)- (b) The linear form i?b'K fc } in 
diagrammatic notation for the single-site case. As before, we leave the bottom tensor uncontracted to 
accentuate the definitions of the C1Z overlaps. 



functional with respect to this block is given by the projected effective Hamiltonian eigenvalue problem 
pbl ttfgpb'Ubl . . . A V+-i] = EA U\ . . . A [j+S ~ 1] , (51) 

where P^ is a projector into the space orthogonal to the {\4>k)}- Given the states as MPSs, 

we construct these projectors as follows. The diagram corresponding to the overlap (ip\4>k) is shown in 
Fig. E(a), with the bold lines corresponding to \<pk) and the thin lines to This is a linear form in 
all of the MPS tensors of and so the condition that be orthogonal to this state for the given 
block of tensors with all others held fixed may thus be stated as 

Mfa) = f ^M* 9 tu^ Mfo)) A® . . . = (52) 



= F iMk} A [j] _ _ _ A b+s-i] = _ ( 53 ^ 

Here the linear form enforcing orthogonality pblw is shown diagrammatically in Fig. Efb). We can 
construct the projector P^ as 

pbl = 1-J2 F [Mk} (-^L, F [Mk ' } f , (54) 
kk' 

where (A/ --1 )^^' is the kk' th element of the inverse of the Gram matrix 

N kk i = Tr (p^« tpblffe'} j . (55) 

This Gram matrix inverse is important to ensure idempotency of the projector. As before, direct 
construction of the projected effective Hamiltonian leads to an algorithm which scales poorly as 0(x 6 )- 
Hence, it is important to use sparse methods which require only the application of P and H on some 
block of tensors A^ . . . A^ +S ~ l \ Direct application of P^ as written requires 0(x i ) operations and also 
scales quadratically in the number of eigenstates desired Ne due to the double sum in Eq. (154j) . To find 
a total of Ne eigenstates by this method thus requires 0(N E x A ) operations, which is unacceptably slow. 
A simple idea to reduce this scaling would be to find the eigenvectors of the inverse Gram matrix and re- 
express the projectors P^lw in terms of them, rendering the double sum a single sum. However, while 
the Gram matrix M is Hermitian and positive semidefinite it may also be badly conditioned and singular. 
A numerically stable alternative to this idea is to construct the Moore- Penrose pseudoinverse [3H] of the 
Gram matrix 

n p 

a=l » 
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where V is the matrix with the eigenvectors of M as columns and A are the n p eigenvalues of M which 
are greater than nA/A max e, where n is the linear dimension of TV, A max its largest eigenvalue, and e the 
machine precision. We use this pseudoinverse to transform to a new set of linear forms 

= J_^^ F b1« (57) 



A^ 



such that 



pm = i-J2 G [Mfl} G mil} f . (5? 



/' 



Often, the dimension of the set {G^^} is much smaller than Ne- The diagonalization of the Gram 
matrix requires 0(N E ) operations, independent of Xi an d its construction and the construction of G in 
Eq. (|57|) both require 0(x 2 ) operations. The operation of on the variational block of tensors is now 

pb]j\b] . , . = _ _ _ j^j+B-i] _ Tr (G [i]{M} f v4 [i] . . . A [j+S ~ 1] ) G [j]{,l} (59) 

which is linear in Ne and scales only as 0(x 2 )- Thus, the dominant scaling for typical parameters 
X S> Ne is still the 0(x 3 ) scaling of the effective Hamiltonian multiply, and the algorithm to find iV^ 
excited states scales as 0(NeX 3 )- A sweeping approach is used as in the ground state search algorithm, 
and the iteration is stopped when the variance drops below a user-specified tolerance. As before, the 
variance does not guarantee that the state found is the next lowest-lying eigenstate, but this does not 
usually cause problems in practice. 

As with the LR overlaps used in the variational ground state search, one can also cache the overlaps 
CTZ used to construct the linear forms F using the recursions 

42 fc} = <w<w , (so) 

£[^+i]W _ A [i]if C [e]k B [i]{k}i (61) 

i 

and 

Tlf^ 1]k = , (62) 

ftMW = B [i]{k}i lZ [i+1]{k} A [e]if (63) 

i 

where BMW j s the MPS tensor of \4>k) at site £. The linear forms are constructed using these overlaps 
as 

F \j]{k}ij-ij+ s -l _ V r\j\{k}ij DblW% R[i+s-l]{k}i j+s -i V \j+s]{k} (aA\ 
r aj3 — 2-^1 "71 -°7i72 • • • -°7s7s+i 7»+i/3 ' ^ > 

7i-7s+i 

see Fig. [5]^b). 

The variational ground state algorithm presented above is essentially equivalent to standard DMRG, 
aside from the calculation of the variance [12]. When finding excited states, however, DMRG-based 
approaches target multiple excited states in a single MPS, which causes the bond dimensions to grow 
and the quality of each individual eigenstate to degrade. Furthermore, as the ground state and all 
excited states are solved together in that approach, the sparse eigensolver must be able to converge 
interior eigenvalues, which is known to be troublesome [381 EH] . We call the present algorithm, which is 
a sparse and numerically stable variant of that proposed in Ref. [19] for PBC, eMPS to accentuate the 
difference. 
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In our experience, there are two main limitations of eMPS. The first is that it is difficult to construct 
good variational guesses for the excited states in contrast to the ground state where the infinite size 
MPS algorithm (251 [50] is applicable. Here, the usefulness of the variance becomes readily apparent, as 
the discarded weight can be 10 -12 or less while the variance is of order 10 -2 in early sweeps. The second 
is that that the area law considerations which make MPS algorithms so practical for ground states do 
not in general apply to bulk eigenstates, and so the bond dimension required to accurately represent a 
general eigenstate may be exponential in the system size, rendering eMPS inapplicable. 

The ability to find excited states is useful in many contexts. It provides access to the dynamical 
gap for determining the location of second order quantum phase transitions [STJ and Kibble-Zurek 
scalings, even when the gap is not between different symmetry sectors. It can help in understanding the 
structure of conformal field theories by providing access to the primary scaling fields [52]. Excited 
states yield the structure function and other dynamic response functions of low-lying excitations. 
Such response functions are of great use for comparing to experimental measurements. Finally, by 
considering more complex functionals such as (ip\(H — e) 2 \ip) — for minimization, one can 

determine level spacing statistics in a desired energy range for systems much larger than are amenable 
to exact diagonalization. Such studies are immensely useful in discussions of integrability and quantum 
chaos, as well as investigations of the thermalization hypothesis [531 EH 155] . 

5. Calculation of observables 

We now turn to how we can extract information from a state expressed as an MPS. We do so by 
the expectation values of Hermitian operators, or observables. We will demonstrate how to compute 
observables of three different types: local observables, two-point correlation functions, and general 
MPOs. 

We define a local observable as an operator which acts only on the Hilbert space of a single site: 
Q[k] _ ^ ii ,of)\i)(i'\. If this site corresponds to the orthogonality center of then the expectation 
value reduces to 



without loss of generality. If the orthogonality center of the MPS, k, lies within the range q < k < r, 
then we can evaluate the expectation value using only the tensors in this range. The most efficient way 
to proceed is first to form the matrix 




iki' h 

The overall scaling for fixed site index k is O (x 2 d 2 ). 

A two-point correlation function is an expectation value of the form (0^0^) where we take q < r 




(66) 



recursively generate 




(67) 



for t = 1, . . . , r — q + 1, and then evaluate 





For fixed q and r, the algorithm scales as O (x 3 ^ + X 2 d 2 )- 
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To compute the expectation of a general many-body observable O expressed as an MPO we start 
from the right (left) boundary and follow the recursion Eq. (I46p (Eq. (148 p ) all the way to the opposite 
boundary, at which point the remaining lxlxl tensor is the expectation value. The overall scaling 
is O (x 3 dxo + X 2 d 2 M (xo)) where xo is the bond dimension of the MPO and M (xo) is the number of 
nonzero elements in the MPO, which usually scales as xo- 

Entanglement measures such as the bond entropy 

Sj = -Tr pj log pj , (69) 
Pj = Tr i<j+1 \iP)(iP\, (70) 

can be calculated from the singular values S of the MPS tensor as 

^' = -E S 7 l0 g S 7 ( 71 ) 

7 

when this tensor is the orthogonality center. These singular values are computed automatically as part 
of the algorithm to shift the orthogonality center, see Sec. 12.11 



6. Time evolution with MPSs 



We now turn our attention to a variational solution of the time-dependent Schrodinger equation 



H(t) \^(t)) 



(72) 



using MPSs. The general strategy is to find some representation of the propagator over some time 
interval [t, t + St], U(t,t + St), and variationally optimize the functional 

2 



\jP(t + 5t))-U(t,t + 5t)\jp(t)) 



(73) 



with respect to the MPS tensors of \tp(t + St)). Several complications arise in this case which were 
not present in the earlier algorithms. The first practical consideration is that the MPO form of the 
propagator may be difficult and very expensive to calculate. The second difficulty is more physical; the 
time-dependent state following a global quench of a Hamiltonian parameter has entanglement which 
generally grows linearly in time [56]. This causes the bond dimension x to grow exponentially in time, 
and so there is some finite time where an MPS simulation will exhaust the available computational 
resources. However, many important questions regarding non-equilibrium dynamics can still be 
answered by considering moderately sized systems and short times. In addition, consideration of a 
situation in which the Hamiltonian changes only locally can greatly increase the accessible system sizes 
and simulation times [53, [58j 159] . 

The most common approach to time evolution for MPSs is to use the Suzuki- Trotter expansion 



cxp 



'L-l 



-iSt H„ 



, n=l 



L-l 

Jexp 



n=l 



2 r 



n 

n=L-l 



exp 



iSt 



O (St 3 ) 



(74) 



or its higher order variants to construct a series of two-site propagators which can be constructed and 
applied easily. This is the basis of the equivalent [T5] TEBD [6U] and tDMRG [SH E2] algorithms. Here 
H n is the nearest-neighbor bond Hamiltonian acting on sites n and n + 1. This approach is no longer 
viable when the Hamiltonian has longer-ranged terms, and attempts to accommodate such longer- 
ranged terms often exhibit poor scaling [S5| I6~4"] and require Hamiltonian-specialized implementation, 
resulting in inefficient, sometimes prohibitively inefficient code. Krylov-based time evolution, which will 
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form the basis for our approach, has been considered in both DMRG [65] and MPS [13] variants for the 
time- independent case. We note that the latter approach has been used [T3], EE] to study time dependent 
systems, but this necessitated very small time steps set by the rate of change of the Hamiltonian in 
order to provide accurate results. Our approach generalizes the latter method to the time-dependent 
case where the error is independent of the rate of change of the Hamiltonian and demonstrates how the 
algorithm can be formulated entirely in terms of FSA rules for MPOs. 



6.1. Commutator- free Magnus Expansions 



The propagator of a general time-dependent Hamiltonian which does not commute with itself at different 
times is given as a time-ordered exponential 

ft+<5* 



U(t,t + St) = T 



cxp 



dt'H (t'] 



(75) 



t+st 



t+st 



whose most well-known form is the Dyson series 

/t+st 
dt 1 H(t 1 ) + J dhj dt 2 H Ox) H (t 2 ) + . . . (76) 

This formulation of the propagator is not convenient numerically, as the Dyson series is an asymptotic 
series and so it can be difficult to determine an appropriate criteria for termination of the series. 
Furthermore, keeping only a finite number of terms in the Dyson series does not preserve the Lie group 
structure of the propagator; that is, the finite approximation is not unitary. An alternative approach 
which produces unitary approximations to the propagator was given by Magnus [67] who used the 
ansatz 



U(0,t) = exp (-it(l (t) 
to define the Magnus series 

oo 

Q(t) = Y,tin (t) 

71=1 



(77) 



(7* 



where the n th term is of order t n in the sense that its power series in t starts with t n . The term Ct n (t) 
involves n nested integrations over n — 1 nested commutators of H (t) at different times. Explicitly, the 
first few terms are: 



6i (t) 



^2 (t) = ~ 



dt x E (tj 



dt\ 



n 3 (t) = - dh 

6 Jo 



dtn 



dtn 



(79) 
(80) 



1-2 



dU 



H{h), H(t 2 ),H(U 



While approximations obtained from truncating the series yield unitary propagators, these expressions 
are still formidable numerically, involving nested commutators and multidimensional integrals. The 
commutators pose a special difficulty for MPOs, as exact multiplication of MPOs involves multiplication 
of the bond dimensions of the MPOs and hence the algorithm scales exponentially in the number of 
terms kept in the series. Optimization algorithms which attempt to variationally shrink the bond 
dimension of an MPO sum or product such as those proposed in Ref. [16] may also be used, but these 
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become numerically unstable for large systems, and when MPOs are subtracted as in commutators large 
cancellations can cause these algorithms to become stuck far from the variational optimum. 

Hence, rather than work directly with the Magnus series, Eq. (178]) . we start from ansatze of the 

form 

s 

U(t,t + 8t) = Y[ exp (-iSttt^ (82) 
i=i 

where each one of the Cli is a linear combination of H at different times in the interval [t, t + 8t] , and 
require that our ansatz matches the Magnus expansion (equivalently, the full propagator) up to order 
8t +1 . We will call such an ansatz a commutator-free Magnus expansion (CFME) [681 EH]. This ansatz 
has a number of features which make it desirable for our purposes. It is exactly unitary and so the norm 
is conserved. Also, provided that we consider the case where only the Hamiltonian parameters change 
in time and the operators are time-independent, the sums of the Hamiltonian at different times can be 
represented exactly as an MPO using the rules of Sec. 12.21 Thus, the need for complex operations with 
MPOs vanishes. Finally, because the ansatz takes into account the time dependence of the Hamiltonian 
explicitly, the time step is not necessarily fixed by the rate of variation of the Hamiltonian, allowing for 
more coarse stepping in time with fixed error. 

Following Ref. [70], the procedure for generating an iV th -order CFME is to expand the function 
H (t) in terms of (shifted) Legendre polynomials P n , 

H n = (2n -l)6t f dxH (t + xSt) P„_i (a?) . (83) 



Jo 

The orthogonality properties of the Legendre polynomials allow the nested integration to be done 
exactly, leaving a series of nested commutators of the H n . Working in a Hall basis [71], this series of 
commutators is matched with the original Magnus expansion to yield the order conditions fi, n such that 

N 

Cu = ^2ft, n H n . (84) 

71=1 

We note that these order conditions are independent of H (t) by construction, and so are set by the 
choice of CFME alone. As we only require the result to be valid to order 8t N+1 , the integration required 
for the coefficients H n may be performed using Gauss-Legendre quadrature of order N/2 + 1. The end 
result of the analysis is that an iV th order expansion with s exponentials may be written as 

U(t + St, t) = e~ isCl1 . . . e- isCl ° (85) 

N/2+1 

fii = Yl 9i,mH (t + x m 8t) (86) 

m=l 

N/2+1 

9i,m = W m S ^2 ( 2n ~ 1) Pn-l (X m ) fi,n (87) 

71=1 

where x m and w m are the points and weights for Gauss-Legendre quadrature [72]. In this work we use 
a fourth order expansion with three exponentials (N = 4, s — 3). The optimal order conditions for this 
case, obtained in Ref. [70], are 



A,! = 11/40, f 1>2 = 20/87, /i, 3 = 7/50, / 2>1 = 9/20 , / 2 , 3 = -7/25 , (88) 
with / s _j + i in = (— l) n+1 /j in . Order conditions for higher order expansions may also be found in Ref. [70] . 
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We now consider that our time-dependent Hamiltonian MPO is constructed from a set of FSA 
rules ^TZ P ^ . . . , O p " j , {h p (t)} , w p ^j j in which the OA and the weights are chosen without loss 
of generality to be time-independent. In this case, the expansion Eq. fl83|) is applied individually to 
each Hamiltonian parameter h p (t), resulting in the parameters {h^}. Now, because of the canonical 
decomposition of Sec. I2.2[ the MPO forms of Cli from Eq. (|86|) at time t can be constructed exactly 
using the FSA rule set ^JZ P f{o pi , . . . , 6 Pn \ , g i)m h p (t + x m 8t) } , w p ^j \. We note that each 

one of these operators has the same bond dimensions as the original Hamiltonian, and the updates of 
operators fij & t each time step can be done in O {Lxo) time, which is essentially negligible. 

The fact that we can construct the fij from the same FSA rules as the Hamiltonian implies that 
our CFME ansatz is equivalent to evolving our system according to piecewise constant Hamiltonians 
of the same form but with differing Hamiltonian parameters. Additionally, as also occurs in high-order 
Suzuki- Trotter expansions, evolution backwards in time may occur. Finally, we note that even terms 
in the Hamiltonian whose parameters do not vary in time have their magnitude altered by Eq. (]86p . as 
Em 9i,m ^ 1 in general. 



6.2. Krylov subspace propagation 

Using the CFME Eq. (j85p we never need to explicitly form an MPO representation of the propagator 
provided we can find an MPS representation of the exponential of an MPO applied to an MPS. We find 
such a representation from minimizing functionals of the form 

2 



exp 



-iStn 



where, importantly, Q has a known MPO representation. We do so by forming a Krylov subspace 
approximation to the exponential [73] in which the Krylov vectors are represented as MPSs. Specifically, 
we do so via the Lanczos algorithm for the matrix exponential, which can be stated as 

Input: Input Cl in MPO form, in MPS form, and a tolerance e for truncating the recursion. 
Initialize: Set /3q = = 1 and |r) = \ip). 

Iterate: For j — 1, 2, . . . until convergence 

Vj) =J r >//Vi- 
r) = Q\vj). 



U 
(ii 
fiii 



(b 
(c 
(d 
(c 
(f 
(g 
(h 



"j = {Vj\r). 



Orthogonalize |r) against \vj) and 

Re-orthogonalize \r) against all v^, k < j if necessary. 

Pj = ( r \ r ) 

Form the matrix exponential of and obtain c^' 

Test for convergence. 



jjU) 



(iv) Finalize: Set \<p) = YH=i c f'\ v i) 

Here T^' is the symmetric tridiagonal matrix with the c^, 1 < i < j on the diagonal and 
1 < « < j — 1 on the superdiagonal. It is important to use a matrix exponentiation method which 
produces a unitary matrix to machine precision in order to not lose the Lie group structure. Because 
of the small linear dimensions of the matrix exponentiation by direct diagonalization is practical. 



An a posteriori estimate for convergence of the Lanczos recursion is that 



2/ViC 



U) 



< e, where e is 



the tolerance [73] • This can be compared with residual estimates in the ordinary Lanczos algorithm for 
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finding eigenvalues. A rigorous bound on the approximation ||^)kryiov — |^)| < 12exp 
can be established [H] when n > p5t/2 with n the number of Lanczos vectors and p 
the spectral width of fl. This estimate shows that for typical tolerances e = 1CT 6 to 10~ 10 , 6 to 20 
Lanczos vectors suffice. 

As stated before, MPSs do not form a vector space and so the multiplication by Cl, the 
orthogonalization, and the final summation cannot be done exactly while keeping the bond dimension 
of our MPS fixed. However, just as with the eigenstate search, we can devise variational algorithms for 
these three operations which are iteratively performed until a desired tolerance is reached and use this 
tolerance to bound the bond dimension of our time-evolved MPS. We begin by briefly reviewing the 
standard algorithm [12] for finding the optimal MPS |0) representing a sum of MPSs X^b=i c fclV'fc) to 
a given tolerance, as the other algorithms are similar but more complex. In this case we have a set of 
CR overlaps defined between our variational state (<f>\ and the states \ip k ) as in Eqs. ( 161]) and ( 163]) . We 
now sweep through the lattice and make the replacement 

^1...#+-i] = ^ Cj , j pb1W ) (90) 



where the F^ k ^ are formed as in Eq. ( 164]) . see also Fig. |5^b). The orthogonality center of |0) is then 
shifted, the CTZ overlaps updated, and sweeping continued. Convergence can be monitored via 



Tr \(A [j] . . . A [j+S ~ 1] y A [j] . . . A [j+S ~ 1] 



(91) 



- 2 ^ Re Tr \c k (A® . . . A^ 8 ' 1 ^ 



< £ 



with Re denoting the real part. Because we do not have to solve an eigenequation at each iteration, 
this algorithm is often much less costly than the iterative eigenstate search. Also, when we have that 
the coefficient vector c and all of the {(^fc)} have length 1, we can normalize the state \4>) at the end of 
the calculation if required. 

The algorithm to variationally fit an MPS to is similar. In this case we have a set of LR 
overlaps defined via the recursions 



-[i] 

' 'recta' 
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(92) 
(93) 

(94) 
(95) 



where the MPS tensors of \<p) are denoted by A and those of \ip) denoted by B. We now sweep through 
the lattice and make the replacement 

A m . . . A \j+s-i] = h^B® . . . B^ +s ~^ (96) 

where the effective Hamiltonian is formed from the LR overlaps as in Eq. ([43]) . Convergence can be 
monitored via 



H\ 



Tr 



(A [j] . . . A [j+S ~ 1] ) j A [j] . . . A [j+ ' 



-i] 



\H 2 



(97) 
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Here, {ip\H 2 \ip) can be computed in a manner similar to the variance, and need only be computed once 
at the beginning of the calculation. We have also assumed that the block of tensors in Eq. (1971) contains 
the orthogonality center. Again, this algorithm is often much less costly than the iterative eigenstate 
search. 

We now turn to steps (iii) (d) and (iii) (e) of our algorithm. Step (iii) (d) is usually stated for ordinary 
vector spaces as 

\r) = \ r ) - aj\vj) - fc-^Vj-x) , (98) 

as a, = (vj\r) and = (vj_ 1 \r) and so this is equivalent to classical Gram-Schmidt orthogonalization. 
Hence, we could implement step (iii) (d) by using the fitting algorithm Eq. (|9~U1) to find the MPS closest 
to |r) — ctj\vj) — /Sj-il^j-i). However, we have found that the following algorithm, which is closely 
related to eMPS, often converges more quickly and also is applicable to step (iii) (e) . In our method 
we look for the optimal MPS \<p) representing but also subject to the constraints that (4>\ipk) — 
for some set {iV'fc)}- We start by copying the state \ip) to a variational ansatz \<f)). We then construct 
overlaps between the state (0| and which we call LI% and a set of overlaps between (0| and \ipk), 
which we call £7Z. We then sweep through the lattice and make the replacement 

a ...a — 7 j iv Q7i d 1i12 . . . £> 7s7s+1 ^t 7s+1/ g , lyjj 

71— 7s+l 

with B the MPS tensors of and A the MPS tensors of \<f>). We then apply the projector into the space 
orthogonal to the \ipk) by constructing the set {Cr^^H, jj, — 1, . . . ,p, as in Eq. (l57j) and performing 

A \j] _ _ _ A v+-i] = A m . . . A \j+s-i] _ Tr (Q[M^ A [j] _ _ _ A b+s-i]^ G ym ^ ( 100 ^ 

for \i — 1, . . . , p. Using the fitting algorithm of Eq. (|9"U|) corresponds to replacing the £71 overlaps, which 
are between the variational state (0| and the set with a set of £7Z overlaps between the state 

and the set Our algorithm, which amounts to fitting followed by modified Gram-Schmidt, uses 

information about the distance between the variational state and those to be orthogonalized against to 
determine operations, and hence often converges more quickly and is more stable. Convergence can be 
monitored by ensuring that (4>\ipk) are orthogonal to a precision e via 

TrU^F^ k A\<e. (101) 

If one requires additional truncation of the bond dimension, one can switch to the ordinary fitting 
algorithm Eq. fl90|) at this point, using a new variational state |£) to fit to \<p). 

We now pause to consider the sources of error in the time-propagation routine. First, because the 
CFME expansion Eq. (|85|) is only of order St N+1 , the error incurred in using this form of the propagator 
is ccfme = ctfi nSL i5t N , where tfi na i is the final time reached. Thus, as the final time desired becomes 
longer, smaller time steps should be taken in order to keep the error fixed. The coefficient c can be 
determined by using this known scaling and decreasing the time step by a constant factor. Factors 
in the range 2 l ' N to 3 X ' N are practical. This strategy can also be used to devise adaptive time-step 
applications such as those used widely in ordinary differential equation solvers. Next, there is the error 
^Kryiov incurred in the Krylov subspace approximation to the exponential. As discussed above, this 
error can be minimized by adding more and more Lanczos vectors. This error also involves the time 
step 5t, and so when determining the CFME expansion error constant c one should be careful that 
fKryiov < ccfme- Finally, there are errors resulting from the variational fitting of MPSs in steps (iii)(b), 

jj In this initialization all of the LR are Kroncckcr deltas provided that \ip) has an orthogonality center. 
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(iii)(d), and (iii) (e) of the Lanczos algorithm for the matrix exponential. These can be reduced by 
lowering the variational tolerances, but this is done typically at the expense of a larger bond dimension 
X and hence a slower algorithm and more memory usage. 



6.3. Simulation Protocol 

We are now in a position to devise a complete, generic protocol for the time evolution of a ID quantum 
system. 

(i) Input: Input an operator alphabet and a set of FSA rules defining the Hamiltonian MPO. Input the 
functional forms {h p (t)} of its Hamiltonian parameters, a final time desired tfi na i and a time-step 
St. Input tolerances {e} for variational ground state search and time evolution. 

(ii) Initial state preparation: Find the ground state of the Hamiltonian using variational ground state 
search from Sec. [3j Alternatively, if a different initial state is desired, read in its MPS representation. 
Set \tp (t = 0)) to be this state. 

(iii) Measure: Measure local observables and two-point correlation functions constructed from the OA 
as well as MPO observables constructed from their own FSA rule sets. 

(iv) eMPS: Use eMPS (Sec. H]) to find excited states of the Hamiltonian, if desired. Measure properties 
of the excited states. 

(v) tMPS: Set t = 0. 

(a) For i — s, . . . 1, construct f2j from Eq. f )86|) using the FSA rules and use the Krylov algorithm 
to apply the matrix exponentials of these operators in succession to \ip (t)). Set t = t + 5t. 

(b) Measure \ip (t)) if desired. 

(c) Continue (a) and (b) until ta na \ is reached. 

Starting from a base set of rules such as the site, bond, exponential, decaying function, and finite 
function rules of Sec. 12. 2\ a single implementation can accommodate a vast range of systems based on 
the particular OA, rules, and quench protocols {h p (t)} used. 



7. Case studies 



7.1. Case Study: Ising model in a transverse field 

For our first case study we choose the paradigmatic Ising model in a transverse field, with Hamiltonian 

H=-jJ2t?ZUi-hJ2*?- ( 102 ) 

Here J is the coupling energy, h is a transverse magnetic field, and the {<jj} are the Pauli spin operators 
on site i. We choose this model because its dynamics are amenable to numerically exact study using the 



time-dependent Bogoliubov-de Gennes formalism (see Appendix A for a review) and so we were able 
to carefully check convergence of our results. The statics have all been converged to eight digits, and 
the dynamics at all times to at least four digits. Here we refer to convergence of local quantities such 
as the energy or density of defects. Nonlocal quantities such as the bond entropy will not have this 
same accuracy, but numerical tests show that the qualitative behavior is unaffected. We begin with a 
discussion of the statics. 

In Fig. [6]^a) we demonstrate the gaps from the ground state to the two lowest eigenstates, computed 
via eMPS. The variances are smaller than the point size in this case. The upper (lower) curve corresponds 
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h/J L 

Figure 6. (color online) Statics of the Ising model, (a) The gaps to the two-lowest-lying eigenstates, 
computed using eMPS. (b) The bond entropy Si vs. the site index i and magnetic field h for L = 100 
sites. Note the increased curvature near the critical point h = 1. (c) The central charge c (blue line) 
and bulk entanglement a (red line) extracted from a fit to Eq. (|103|) vs. h for L = 100 sites, neglecting 
30 sites at both boundaries. The inset is a close-up of the critical region, (d) The bond entropy of the 
central site Sx/2 vs. L at h = 1.0 together with a fit extracting the central charge. 



to even (odd) parity, while the ground state has even parity, where parity is defined as simultaneous 
inversion of all spins P = (J~J f af). Hence, the relevant gap for discussing the quantum phase transition 
is in fact the gap to the second excited state, shown in greenjjj] Both gaps close at criticality, and this 
can cause the first excited state returned by eMPS to be a mixture of these two states. This will not 
affect the energies so long as the variance remains small, but it can affect other observable properties of 
the states. There are two ways to remove this nuisance. The first, most complex, and most preferable 
is to use an MPS representation in which the state is explicitly Z2 invariant [39]. The second is to add 
a field coupling to the parity — h p Y[ { a* to cause the relevant even-parity subspace to become lower 
in energy than the odd-parity subspace. This operator is an MPO with bond dimension 1, the MPO 
equivalent of a product state. The closing of the gap at the known critical point h — 1 is linear in 1/L, 
indicating a conformally invariant critical theory with dynamical critical exponent z — 1. 

We venture to determine the central charge of the critical theory by fitting to the Calabrese-Cardy 
formula in two ways. In the first, we fit to the finite-size formula 



S< = ^ log 
o 



"2L 


^7rA 


— sin 




7T 





(103) 



ff The even parity gap for ft, > 1 is in fact twice the demonstrated gap in green, but the essential piece is the closing of 
the gap at criticality. 
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for fixed L and variable i, and in the second we fix % at L/2 and fit Sl/2 to this formula for various 
L. Near criticality, the presence of nonzero c indicates a curvature of Si, while in the gapped phases 
Si obeys an area law and is hence a constant apart from finite-size effects. The bond entropy in the 
bulk approaches the correct limiting values of log 2 as h — > and as h — > oo. The first fit, shown in 
Fig. 13(a), provides us with a clear indicator of the critical region by the spike in the central charge c. 
However, the precise determination of c for a finite size system in this case is noisy, likely due to strong 
finite-size effects. Once we have narrowed down where the critical region is, the second fit, shown in 
Fig. 16(b), allows us to extract the anticipated value c = 1/2 much more precisely. If the same scaling 
analysis is attempted at a point which is not the critical point, the bond entropy saturates and c — > 
as L — > oo. We can understand this as a large but finite correlation length £. For L/2 < £, the system 
appears to be conformal and we see scaling of the bond entropy with L. For L > £ the bond entropy 
saturates and this scaling breaks down, indicating that the given region is not critical. We note that in 
this analysis we have used no properties which are specific to this system e.g. correlation functions of 
an order parameter to extract the critical behavior. 

We now turn to the dynamics. The Ising model has also been a subject of interest for dynamics as 
a testbed for the Kibble-Zurek hypothesis (KZH) that equilibrium properties determine nonequilibrium 
properties following a quench across a quantum critical point. This was studied in Refs. [75l 1761 177] 
using the quench protocol 

h (t) I J = 5 (l - , < t < r . (104) 

A useful quantity for determining how non-adiabatic a particular quench is in this case is the density 
of defects 

^^E^-^+O ( 105 ) 

i 

which is the density of Bogoliubov quasiparticles at zero magnetic field. In addition to the density of 
defects, universal scaling has also been predicted in the heat, or non-adiabatic part of the energy, 

Q (t) = (t) \H (t) \i, (t)> - (V g ., (t) \H (t) |V g ., 0)} (106) 

where |^> g . s .(i)) represents the instantaneous ground state of H(t). In addition to these quantities, which 
are amenable to Bogoliubov-de Gennes analysis, we also compute the time-dependent bond entropy. 
Our results are shown in Fig. [7J 

We first discuss the heat, as shown in panels (a) and (d). Panel (a) displays the heat as a function 
of time, and demonstrates a sharp change in the behavior of the system as we pass through the critical 
point. This is especially true of the longest quenches. In panel (d) we investigate the heat as a function 
of the quench rate both at the time when h takes on its critical value, t c = 4r/5, and at the final 
time when h = 0. The large difference indicates that non-adiabatic processes continue after we have 
passed from the critical region back into the gapped ferromagnetic region. Thus, the universal scaling 
of the heat may be difficult to determine if the critical point itself is not known sufficiently well. We 
now move on to the density of defects as shown in panels (b) and (e). In panel (b) we see that the 
density of defects at the final time decreases slowly to zero as r — oo; that is, when the quench becomes 
perfectly adiabatic. This is in accordance with the KZH prediction. Panel (e) demonstrates the large 
disparity between the density of defects at the critical time and the final time for all but the most rapid 
of quenches. Finally, in panels (c) and (f) we investigate the bond entropy. In panel (c) we see the bond 
entropy of the central bond as a function of time. As the quench becomes more adiabatic, the entropy 
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increases more towards the ferromagnetic limiting value of log 2. However, for very slow quenches, the 
bond entropy reaches this value before the end of the quench and then begins to oscillate. In panel (f) 
we show the bond entropy as a function of the bond index at the critical time. Bulk curvature such 
as that seen at criticality in Fig. El^d) is absent, indicating that the time-evolved state is still quite far 
from the conformal ground state. 

1.2. Case Study: Dipolar Ising chain 

In this section we go beyond exactly solvable models and study a dipolar Ising chain 



Such models are relevant to the study of ultracold molecules in optical lattices, where the dipole-dipole 
interaction falls away as 1/r 3 with r the distance between dipoles [TSJ [TJ, [HI El HZ]- Here the cutoff 
\j — i\ < 6 represents a consistent order of approximation in going from a Hubbard- type model with 
dipolar interactions to a spin model. We stress that all results obtained in this section were obtained 
using the same implementation as the last section. 

We first turn to the statics of this model, shown in Fig. [SJ Many of the features are similar to 
those of the nearest-neighbor Ising model. The most important differences are that the critical region 
is shifted towards larger h with respect to the nearest-neighbor Ising model as seen in panels b) and c). 
This indicates increased stability of the ferromagnetic phase, in accordance with expectations. Using 
these points as a guide, we determine the critical point to be h c = 1.362 ± 0.01, as shown by the scaling 
in panel d). Additionally, the energy of the first even parity excited state deviates from pure z = 1 
linear behavior near h = 0, indicating interactions between quasiparticles which were noninteracting in 
the nearest-neighbor Ising model. 

We now turn to dynamics, following the same quenching protocol Eq. (11041) as above. The results 
are shown in Fig. [9j We reiterate that the dynamics of this model cannot be handled by Bogoliubov-de 
Gennes methods, nor straightforwardly with standard tDMRG/TEBD. The density of defects no longer 
represents the density of quasiparticles at the critical point, but we compute it for comparison with 
the results of the nearest-neighbor Ising model. Because of the larger MPO bond dimensions and the 
more rapid growth of bond dimension for this model, we restrict our analysis to short times Jr/h < 5, 
though an optimized implementation could reach longer times. The basic features are similar to the 
dynamics of the nearest-neighbor Ising model. One quantitative difference is that, because the critical 
point is reached at an earlier time than in the nearest-neighbor Ising model, oscillations in the central 
bond entropy occur for shorter quench times. 

8. Conclusions 

The power of matrix product state algorithms over DMRG-based algorithms is most readily apparent 
when multiple states are involved, as each state may be represented as a separate matrix product state 
in the former approach. Because matrix product states with a fixed bond dimension do not form a 
vector space, a set of matrix product states carries more information at smaller numerical cost than 
the same set represented as a multi-state targeted basis in DMRG. We have presented two algorithms, 
eMPS and tMPS, which use this property to find eigenstates and perform time evolution of strongly 
correlated ID quantum systems. 




(107) 
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Figure 7. (color online) Dynamics of the Ising model, (a) The heat as a function of time shows a 
marked change in behavior as we transition past the critical time t/r = 4/5. The numbers indicate the 
value of tJ/H. tJ/H = 0.1 and 0.01 are indistinguishable on the scale of this plot. Both may be considered 
to be instantaneous, (b) The density of defects as a function of time scales to zero as tJ/H — > oo, in 
accordance with the KZH. (c) The bond entropy of the central bond as a function of time approaches 
the limiting value log 2 as the quench becomes more adiabatic. For nearly adiabatic quenches, the bond 
entropy oscillates after the critical point, (d) Scaling of the heat in the final and critical stages with 
the inverse quench time shows marked non-adiabatic processes occurring after passing the critical point, 
(e) Scaling of the density of defects in the final and critical stages with the inverse quench time shows 
non-adiabatic processes after passing the critical point only for slow quenches, (f ) Snapshots of the bond 
entropy at the critical time demonstrate that the system is not generally close to its conformal ground 
state. 



eMPS uses a set of eigenstates stored as separate matrix product states to define a projector 
into the space orthogonal to this set. We use this projector to explicitly orthogonalize a variational 
state against previously determined eigenvectors in order to find excited states. The explicit matrix 
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Figure 8. (color online) Statics of the dipolar Ising model, (a) The gaps to the two-lowest-lying 
eigenstates, computed using eMPS. Here we see a breakdown of the linear dispersion at small h, indicating 
interactions between quasiparticles. (b) The bond entropy Si vs. the site index i and magnetic field h for 
L = 80 sites. The point of greatest curvature is shifted towards larger h with respect to the Ising model. 

(c) The central charge c (blue line) and bulk entanglement a (red line) extracted from a fit to Eq. (|103[) 
vs. h for L = 80 sites, neglecting 30 sites at both boundaries. The inset is a close-up of the critical region. 

(d) The bond entropy of the central site £l/2 vs. L at h = 1.362 together with a fit extracting the central 
charge. 

product state representation allows us to store the excited states much more accurately than with 
standard DMRG, and allows allows us to ensure global orthogonality between the eigenstates to a desired 
tolerance. The variance, which is computed exactly using the matrix product operator representation 
of the Hamiltonian, gives strict error bars on the energies obtained with this procedure. 

tMPS avoids the need for an explicit representation of the propagator by using a commutator- 
free Magnus expansion and then building successive Krylov subspace approximations to the matrix 
exponentials which appear in the expansion. Each vector in the Krylov subspace is stored as a separate 
matrix product state to maximize efficiency Furthermore, the operators f2j have exact representations as 
matrix product operators with the same bond dimension as the Hamiltonian. Our algorithm eliminates 
the need for Hamiltonian-specialized implementation of dynamics. Additionally, by carefully accounting 
for the time dependence of the Hamiltonian with a commutator-free Magnus expansion, the error in our 
algorithm depends only on commutators of the Hamiltonian with itself at different times and not on its 
derivatives. As with eMPS, the errors are rigorously accounted for by considering distance functionals 
with the variational state. 

The matrix product operator forms of ID Hamiltonians can be obtained using a small set of 
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Figure 9. (color online) Dynamics of the dipolar Ising model, (a) The heat as a function of time displays 
a slower buildup in post-critical non-adiabatic effects for longer quenches, as in the nearest-neighbor case, 
(b) The density of domain walls as a function of time is comparable to that for the nearest-neighbor case, 
but no longer has the same interpretation in terms of quasiparticles. (c) The bond entropy of the central 
bond oscillates for the shorter quench time Jr/h = 5 due to the quench passing the critical point sooner 
than in the nearest-neighbor case. 



finite state automaton rules such as the site, bond, and finite function rules. Using matrix product 
operator arithmetic, we can add together the various terms in a Hamiltonian from these rules to form 
a complete canonical MPO representation. This representation allows for templating of Hamiltonians 
which depend only on the type of interactions and not on the microscopic constituents of the lattice 
model. Furthermore, given the time- dependent form of the Hamiltonian parameters, one can use the 
same template to form the operators f2j which appear in tMPS at negligible numerical cost. 

We used our algorithms to study both the nearest-neighbor Ising model in a transverse field and a 
long-range Ising model in a transverse field. By the closing of the gap obtained with eMPS we determined 
that the critical points of these models were conformal, and so we used the Calabrese-Cardy formula for 
the bond entropy of conformal systems to locate the critical point and its associated central charge. The 
known result h c — 1 was verified for the nearest-neighbor case, and the critical point was shifted deeper 
into the paramagnetic region h c = 1.362 ± 0.01 for the long-range case, indicating a stabilization of the 
ferromagnetic phase. We used tMPS to study the dynamics of these models following a linear quench 
of the transverse field from the paramagnetic phase through the critical point into the ferromagnetic 
phase. We saw strong non-adiabatic effects in the heat as quenching continued into the ferromagnetic 
region, scaling of the density of defects consistent with the Kibble-Zurek hypothesis, and the oscillation 
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of the bond entropy near its limiting ferromagnetic phase value for nearly adiabatic quenches. 
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Appendix A. Time evolution of exact solution of ID transverse-field quantum Ising 
model for comparison with tMPS 

The solution of the statics of the transverse-field quantum Ising model is covered in standard texts 
However, for comparison with tMPS we require a description of the dynamics, and so we present the 
dynamical case here. To find the exact solution of the Ising model, 

L 

we affect the Jordan- Wigner transformation 

af = l- Tc\h , (A.2) 

at = -(ci + 4) II i 1 ~ 2£ fe) > ( A - 3 ) 

j<i 

where the fermionic operators q satisfy the anticommutation relations {q, fit} = Sij, {q, cj} = {cj, cj} = 
0. This transforms the Ising model into the fermion Hamiltonian 

L-l L 

k = - J E ( £ i^i + + h - c -) + 2/1 E d ^ - Lh ■ ( A - 4 ) 

8=1 1=1 

As this is a quadratic form in fermion operators, it may be diagonalized by a canonical (Bogoliubov) 
transformation [791 EH] 

fe=i 

which provides the set of Bogoliubov-de Gennes equations 

e k u k = Au k + Bv k (A. 6) 

efcV fc = -Bu k - Aw k (A. 7) 

where u k are the elements of {ui k , i — 1, . . . , L} arranged as a vector and likewise for v k . The matrices 
A and B are real and tridiagonal, with the nonzero matrix elements A^j = 2h, Aij = — J, \i — j\ = 1 
and B iji+ i = — B i+ i ti = — J. The pairs (u ik ,v ik ) with positive energy e k , e\ < e 2 < . . . < e^, and the 
normalization |u^| 2 + |v^| 2 = 1 define quasiparticle operators 

% = u* ik di + v ik c\ (A. 8) 

which bring the Hamiltonian into the diagonal form 

- / 1\ 

^ = E e M^ fc ~2 • (A - 9) 

k=i \ ' 

Corresponding to every such pair is another pair {u ik ,v ik ) = (v ik ,u ik ) with e k = —e k which defines the 
conjugate operator Writing u^ and v k together as a composite vector, the Bogoliubov-de Gennes 
equations take the form of a real symmetric eigenproblem of dimension 2L: 

which can be readily solved using standard eigenvalue methods [8Tj . 
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Figure Al. Demonstration of Bogoliubov-de Gennes method convergence. The errors in the energy 
computed with a given time step St, E& t , are shown as a function of time. The ratio of the errors as 
a function of time is roug hly (0.01/0. 1) 4 = 0.0001, as should be expected for our fourth-order CFME. 
Errors in the density of defects behave similarly. 



Evolution under the fermion Hamiltonian Eq. ( 1A.4P does not preserve the number of fermions Np 
but it does preserve the fermionic parity (—1) Nf . Because the ground state is the Bogoliubov vacuum 
it contains no fermions, and so the first accessible excited state consists of two Bogoliubov excitations, 
one in each of the lowest two modes. The gap between the ground and first excited states is thus ei + e 2 . 

We now consider the Heisenberg equations of motion for the fermi operators 

ih ( tEM = A(t)c(t) + B(t)c f (t), (A.ll) 

(JjV 

where A(t) and B(t) are the time-dependent generalizations of A and B above and c(t) and &(t) are 
the elements of {ci{t)} and {c\(t)}, respectively, arranged as vectors. Because this equation is linear in 
the Fermi operators it may be diagonalized with a time-dependent Bogoliubov transformation 

L 

Ci (t) = ( u i* 7fc + < (<) 7fe) (A.12) 
fe=i 

where Uj(t) and Vj(t) subject to the time-dependent Bogoliubov-de Gennes equations 






H(t) fc , (A. 13) 



A (t) B (t) 
-B(t) -A(t) 

and 7fc and tf k diagonalize the Hamiltonian at the initial time. Equivalently, u and v define time- 
dependent Bogoliubov operators 

Ik (t) = u* k (t) di + i>ik (t) c] (A. 14) 

such that the time-evolved state \ip{t)) is the Bogoliubov vacuum of this set, i.e. l Jk{t)\i } {t)) = 0. To 
compare with the MPS simulations, we note that the energy at time t is 

-^E c *w = 4^ ( ^ ( * )v * ( * ))w( * ) (^(t)) • (A - 15) 

Similarly, the density of defects is 

? w = 2i £<^ {t) 1 1 1 " ^ {t)) ' (A - 16) 

1=1 
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L-i r l 



i=i L fc=i 



«i,fc (*)) «+i,fc 



(*)) 



(A.17) 



Time evolution thus reduces to the solution of a 2L x 2L time-dependent matrix differential equation 
which we solve using a CFME as in Sec. 16.11 Because the dimensions of the system are much smaller 
than those of the associated MPS problem we are able to take very small time steps, and so the results 
obtained in the method may be considered to be numerically exact, see Fig JAll Here we demonstrate 
that the error incurred in the energy as a function of time scales with the infinitesimal time step St as 
St 4 using our fourth-order CMFE. Hence, by decreasing St, any desired degree of accuracy may be met. 
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